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Abstract. The structure of the reconstruction algorithm OPED permits a 
natural way to generate additional data, while still preserving the essential 
feature of the algorithm. This provides a method for image reconstruction 
for limited angel problems. In stead of completing the set of data, the set of 
discrete sine transforms of the data is completed. This is achieved by solving 
systems of linear equations that have, upon choosing appropriate parameters, 
positive definite coefficient matrices. Numerical examples are presented. 



1. Introduction 

Image reconstruction from x-ray data is the central problem of computed tomog- 
raphy (CT). An x-ray data is described by a line integral, called Radon transform, 
of the function that represents the image. A Radon transform of a function / is de- 
noted Tif{9, t) where 9 and t are parameters in the line equation cos 6'x-l-sin 9y = t. 
The image reconstruction means to recover the function from a set of line integrals 
by an approximation procedure, the reconstruction algorithm. For further back- 
ground we refer to [SJ [SI US]- The quality of the reconstruction depends on how 
much x-ray data is available and the data geometry, meaning the distribution of 
the available x-ray lines, as well as on the algorithm being used. The ideal case 
is when the available data are exactly what the reconstruction algorithm need. 
Most of the algorithms, for example the FBP (filtered backprojection) algorithm, 
requires a full set of data that are well distributed in directions along a full circle of 
views. In many practical cases, however, x-rays in some of the directions could be 
missing. We then face the problem of reconstructing an image from a set of incom- 
plete data, which is, however, intrinsically ill-posed. In order to apply an algorithm 
that requires a full set of data on the problem of incomplete data, one needs to 
derive approximations of the missing data from the available data, for example, by 
some type of interpolation process, which, however, has to be done carefully as the 
incomplete data is usually severely ill-posed. 

In the present paper we consider the limited angle problem, a type of incomplete 
data problem for which the radon data TZf{9, t) are given for in a subset of a half 
circle, and show that the reconstruction algorithm OPED (based on Orthogonal 
Polynomial Expansion on the Disk), studied recently in [T9j [2Qj [21], permits a 
natural approximation for the missing data. The limited angle problem was studied 
extensively in [2J [5J [HI UHl [HI [IS] , see also [13] ■ The problem is known to be highly 
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ill-posed ( [2] ) . The approach in jHl |9j [TOj [TT] uses the singular value decomposition 
to generate the missing data, then uses FBP to reconstruct the image. 

In our approach, we do not actually generate the missing Radon data per se, but 
what is missing for the OPED algorithm, which are the discrete sine transforms of 
the missing data. This algorithm for two dimensional images is based on orthogonal 
expansion on the disk; in fact, it is a discretization of the A'^-th partial sum of the 
Fourier expansion in orthogonal polynomials on the disk. One of the essential 
features of the algorithm is its preservation of polynomials of high degree. In other 
words, if the function that represents an image happens to be a polynomial of degree 
no more than N, then the algorithm reproduces the image exactly. For smooth 
functions, this ensures that OPED algorithm has a high order of convergence. In 
fact it is proved in fl9J that it converges uniformly on the unit disk for functions 
that has second order continuous derivatives. Furthermore, numerical tests have 
shown that the algorithm reconstructs images accurately with high resolution for 
both phantom data and real data. Our main result in Section 3 shows that we 
can make use of the structure of the approximating function in OPED algorithm 
to generate what is missing for the algorithm, while still maintaining the feature 
of polynomial preserving, so that the algorithm can be used for the limited angle 
problem. The method completes the set of discrete sine transforms of the data by 
solving systems linear equations. We show how to choose parameters so that these 
matrices are positive definite. The ill-posedness of the limited angle problem is 
reflected in the ill-conditioning of the matrices. We discuss the dependence of the 
condition numbers on the parameters that appear in the algorithm, which serves 
as a guidance for the numerical experiments. 

The paper is organized as follows. The follows section contains the background on 
OPED algorithm. In Section 3, we derive the algorithm for limited angle problem, 
provide a theoretic background, discuss conditions for the matrices to be positive 
definite, and study the conditional numbers of the matrices. The numerical results 
are reported and discussed in Section 4. A shot conclusion finishes the paper in 
Section 5. 

2. Background and OPED algorithm 

2.1. Background. Let f{x, y) be a function defined on the unit disk B = {{x, y) : 
+ 1^ i}- A Radon transform of / is a line integral. 



where l{0,t) — {{x,y) : xcosO + ysm9 — <} n i? is a line segment inside B. 
The central problem in CT is to recover the function f(x,y), which represents an 
image, from its Radon transforms, which represent x-rays in mathematical terms. 
In reality, only a finite collection of x-ray data is available for reconstruction, which 
can be used to construct, in general, an approximation of /. An algorithm is a 
specific approximation process to / based on the finite collection of data. There 
are many ways to construct the approximation process. The FBP algorithm is 
based on an interaction between Fourier and Radon transforms. OPED algorithm 
is based on orthogonal expansion on the disk. 

Let denote the space of polynomials of total degree at most n in two variables. 
Let Vn{B) denote the space of orthogonal polynomials of degree n on B with 




< 6* < 27r, -1 < t < 1, 
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respect to the Lebesgue measure. A function in L^{B) can be expanded in terms 
of orthogonal polynomials, that is, 

oo 

(2.1) fix) = ^proj, fix), proj, : L^iB) ^ V„(S). 

fe=0 

It turns out that the projection operator projj, / has a natural connection to the 
Radon transforms. In fact, the following expression holds (|I2], see also [71 [HI \T\ ) . 

1 1 

(2.2) projfc /(x, y) - - V - / 7^/(</),, i)t/fc(^)di(fc + l)[/fc(xcos + y sin 0,), 



where (/)^ — and Ukit) denotes the Chebyshev polynomial of the second kind, 

/ ^ sin(fc + 1)6* 
2.3 Ukit) = \ „ , t = cose. 

sm a 



The formula ( 2.2 1 allows us to construct a number of approximation processes based 
on the Radon data. Here are two that are of particular interests to us, 

Af-l N-l 

(2.4) Snfix) ^ proj^ /(x, y) and S^fix) ^ ?7(^) Projfc /(^^^ 2/)> 

fc=0 fc=0 

where 77 is a smooth function in C^[0,cx)) such that rjit) = 1 for t e [0,t], where r 
is fixed with < r < 1, 'qit) = for t > 1, and -qit) is strictly decreasing on [r, 1]. 
The function Sat/ is the best approximation to / from 11^ in L^iB) and it is a 
projection operator on 11^, that is, Sat/ = / if / e 11^, while the function Sjff 
approximates / in uniform norm with the error of approximation in proportion 
to the best uniform approximation by polynomials of degree \tN\ and it satisfies 
S'lif = f ii f e ^IrN] (s*3e [IH])- We can discretize S'at/ or S'l^f, by applying 
a quadrature formula on the integral over t in (2.2 1, to get an approximation to 
/ based on discrete Radon data, which is the essence of the OPED algorithm. If 
we choose Gaussian quadrature with respect to the Chebyshev weight, then the 
discretized approximation functions, denoted by Aj^f or Aj^f, respectively, also 
preserve polynomials of appropriate degrees. 

To be more precise, we work with the following explicit OPED algorithm. 

Algorithm 2.1. OPED Algorithm. Let Nd and N be two positive integers and 
N(i < N . Evaluate at each reconstruction points, 

I N,-1N-1 . ^ . 

(2.5) yljv(a;,y) = — ^ ^ [ — J Afc,,,(/c + l);7fe(xcoS(?!)^ + ysin(/)^) 
where (f)y — 

1 "^^F' • , -r,.. , ^ , (2j + 1)^ 



k=0 1^=0 ^ ° 



(2.6) Afc,^ = — ^ sin(fc+ l)-0J7^((/)l,,cos■0j), ipj 

d A n 



3=0 

and rjit) is a smooth function such that rjit) = I on [0, r] for a fixed r, < t < 1, 
and rjit) > for t > t. 
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The image is reconstructed by the values of ANix, y) over a grid of reconstruction 
points. The function AN{x,y) is a polynomial of degree Nj^. As an operator, it 
preserves polynomials of degree [riV^J, that is, 

ANl = f forall/en^^^^j. 

Naturally Nd and N could be the same. For image reconstruction, we often take N 
and Nd as large as 1000, meaning that Anf preserves polynomials of high degrees. 
The reconstruction has high quality, as supported by both theoretic study in [19] 
and by numerical experiments in [31[2ni[H]- A fast implementation of the algorithm 
is discussed in [20 , which shows that we need 0{N^) evaluations for reconstructing 
an image on a M x M grid, if Nd k, M k, N . 

2.2. OPED algorithm with odd number of views. An x-ray enters an area 
in the angle (j) is the same as the x-ray that exits with the angle vr + 0. For Radon 
transform, this is stated as 

(2.7) 7^(0^- 7r,t) =7^(0,-^), O<0<27r. 

As a result, we have been using the OPED algorithm with N being an odd integer 
to avoid the repetition. For N being odd, we can rewrite the formula of OPED 
algorithm so that the views are restricted to [0,7r] instead of [0,27r]. We state this 
as a proposition. 

Proposition 2.2. Let N be an odd integer. Then we can replace (pu — ^ixv jN in 
( [2l5| and ([2^ by 7^ = TxvjN. 

Proof. Let us define 

^fe(0) = ^ X! sin(fc-|- l)V'j7e(,/,,cosV'j). 

Then X^^u = ^k{4'y)- Since N is an odd integer, it follows readily that satisfies 
4>v+{N+i)/2 — + l2v+i- We also have that ipj satisfies tt — ipj = ipNd-j-i- ^ 
result, it follows from (2.7) that 

T^{(t)u+{N+l)/2,C0Stpj) = 7^(72l.+l,-COsV'j) = 7^(72i.+i,COS?/'Ar^_j_i). 

Then, for <v <{N ~3)/2, we obtain 

Afe(0i.+Ar/2) = ^ sm{k + l)ipjTZ{'j2u+i,cosi;Ni-j-i) 

3=0 



— sin(/c-|- l)V'Ar<j-j-i7?.(72i/-Hi,cosV'j) 



= (-1)'']^ J2 sin(fc -Hl)^J7^(72.+l, cos V,) = (-!)"■ Afe(72.+i). 

j=0 

Let ilki4>) ■= Xk{(t')Uk{x cos (/) + y sin (j)). Using cos cj)^, ^ w+i = — cos72;y+i and 
sin f/)^. I w+i ~ — sin 72^+1, as well as Uk{—t) = {—\)^Uk{t), it follows that 

^k{'l>u+!:^+^) = (-l)''Afe(72^+i)C/fe(-.TCOS72i.+i - ?;sin72i.+i) = VLk{-f2u+i)- 
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Consequently, we obtain 

N-l JV-1 

Xk,vUk{x cos + ysin(j) 



JV-l 



i/=0 i/=0 i/=0 

from which the proof of the stated result follows immediately. 



□ 



2.3. OPED algorithm with even number of views. If N is even, the relation 
(2.7 1 shows that some of the rays coincide, so that the formulas in the OPED 
algorithm can be simplified somewhat. We summarize the essential part in the 
following proposition. 



Proposition 2.3. Ltt N be an even integer. Then Xk,u defined in (2.6 1 satisfy 
(2.8) \k,.+N/2 = (-l)'Afc,,, 0<iy< N/2-1 

and, furthermore, 



N-l N/2-1 

(2.9) ^k,Mx cos(/)y + ?/sin0^) = 2 5^ Xk^uUk{x cos ■ 



u=0 



v=0 



Proof. Since N is an even integer, (p^ satisfies (j)u+N/2 = 
TT — tpj = ^/jTVd-j-i- As a result, it follows from (2.7 1 that 



2/ sin 0,,). 
We still have 



T^{4'u+N/2,cos^j) = 7?.((/)^, -cost/ij) = 7^((/)^,cosV'^f<^-J-l)• 
Following the same line of the proof in the previous proposition, the above rela- 
tion leads to (2.8 1 and (2.9 1 Following the same line of the proof in the previ- 

Similarly, we have in this case 



ous proposition, the above relation leads to (2. 
^ki'Pi^+n) = ^ki'Pi'), from which (2.9l follows. 



As a of consequence of this proposition, the algorithm for even N becomes: 



□ 



Algorithm 2.4. (OPED Algorithm for even N). Let N be an even integer. 
Evaluate at each reconstruction points. 



(2.10) AN{x,y) 



N 



Na-lN/2-l 

E E ' 

k=0 v=0 



k 



Xk,v{k + l)Uk{xcos 4>i, + ?/sin(/i^). 



where (f>i, — and Xk,v o,re given in (2.6 I. 



In other words, we have (2.5 1 replaced by (2.101. Notice that the view angles 
in (2.10) are equally distributed over an half circle; that is, (f)^ in (2.101 are in 
[0,7r]. When we work with the limited angle problem, we will further assume that 
Nd = N/2 in ( |2.10[ ); see Section 4. 

For N being even, a full data set for the OPED algorithm is then 

(2.11) Vn {gy,j := TZ{(f>^, cos tfjj) : < < N/2 - 1, < j < iV - 1} , 
with angle (f>^ distributed equally over a half circle (an arc of 180°). 
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3. Derivation of OPED algorithm for limited angle problem 

In the hmited angle problem, the data available consists of g^j with (j)^ dis- 
tributed over an arc of less than 180°. We are particularly interested in the case 
that N is even and the data is given by 

(3.1) Vr.N ■■= {gu,j ■■r<iy<N/2~l,0<j<N~l}, 

where r is a positive integer and r < N/2—1. In other words, the Radon projections 
correspond to the angles (f'uo, ■ ■ ■ i 4>ur-i ^^'^ missing from the data set I?jv In this 
section we show how the structure of Aj^f can be explored to deal with such a 
problem. 

3.1. Description of the idea. From the given data, we can compute (via FFT) 
every element in the set 

(3.2) A^,A, {Afe,^ : r < 1/ < N/2 - I, < k < Nd - 1} . 

To apply OPED algorithm, the missing data Xk.u for < fc < A^^; — 1 and < ly < 
r — 1 are needed. We now describe our approach to complete the data set. 

Note that the evaluation of AN{x,y) in (2.101 can be carried out so long as we 
know all Xk,^ for < < N/2 - 1 and < fc < iV^ - 1. The equation ( |2.10| is 
derived from (2.5 1 when N is even. For more generality, we work in the following 
with (2.5 1 in which N can be either even or odd, and accordingly with the available 
Xk,i, given by 

(3.3) Ar^N ■■= {Afe,^ : r<i'<N-l,0<k<Nd-l}. 

We will need a lemma on the Radon transform of orthogonal polynomials. 

Lemma 3.1. [12 If P is an orthogonal polynomial in Vk{B), then for each t € 
(-1, 1) andO <e < 2tt, 

nP{9,t) = -^^/l^Uk{t)P{cosB,sme). 

K ^ 1 



Our new algorithm is based on following observation on Xk.u defined in (2.6 1. 

Proposition 3.2. // / is a polynomial of degree at most tN < Nd, then Xk.v 
defined in (2.6) satisfies the system of equations 

forO<k<Nd-l andO < fj. < N - I. 

Proof. If / is a polynomial of degree < tN, then Af = f and we have 

^ Na-lN-l / k \ 

(3.4) f{x,y) = Af{x,y) = — ^ ^ Afe^^r; f — j (fc + l)J7A;(.TCOs</«^ + ysin0^). 

Since TZf{(j}, t) /\J\ — t'^ is a polynomial in t of degree at most riV, as can be seen 
from Lemma |3.1| and we derived (2.5 1 by applying Gaussian quadrature of degree 
2Nd — 1 with respect to the Chebyshev weight, it follows that 

^fe-A- = ixT I] sin(fc + l)?^J7^/(0.,cos?^,) = - / nf{c^,,t)Uk{t)dt. 
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It is known that Uk{x cos (j)^+y sin is an orthogonal polynomial in Vk{B). Hence, 
applying Radon transform on ( |3.4[ ) and using Lemma 3.1 we obtain that 

Integrating against Uk{s)ds and using the orthogonality of Uk, we end up with 



7^/(0,s)C/fc(s)ds = r, 



k 



1 

TV ^ 

1^=0 



Afe,,.?7fc(cos(0 - 



Setting (j) — in the above relation proves the stated relation. 



□ 



Assuming that we are given the incomplete data (3.11. Then we can compute 
Xk^fj, in Ar_N defined in (3.3 1. In order to apply the OPED algorithm, we do not 
need to know each individual missing data. It is sufficient to find the missing Xk^i^', 
that is, to find 

{Afe^^ : < J/ < r - 1, < fc < iVrf - 1}. 
The proposition suggests that we solve these A^^^ from the following linear system 
of equations: For fc = 0, 1, . . . , Nd — 1, solve 



(3.5) 



where for k 



4i = ^ 



Afc,M 

0,1, 
k 



i/=0 



(k) X 



< /i < r - 1, 



. . . , Nd — 1 and < i^, /.i < N ~ 1, we define 
sin(fc+ 1)(0^, - (f>^) 
A^sin(0p - (f)^) 



v ^ fi, and 



k+1 
N 



Notice that A^^^ in the right hand side of (3.5) can be computed from the data in 
(3.11 by (2.6 1, so that they are known. 



To summarize, the idea for the new algorithm is to solve (3.5 I for the missing 
\k,vj OLi^d then apply OPED algorithm to the full set of Xi^ ,^ for reconstruction. 



Solving (3.51 amounts to solve Nd linear systems of equations of size r x r. In 
order for this proposed method to work, it is necessary that the coefficient matrices 
of these systems are invertible, which we study in the following subsection. 

3.2. Non-singularity of the matrices. In this section we assume Nd — N. We 
consider the case that ri{t) = 1 first and define 

, _ sin(fc + 1)(0^ - 0^) 



B 



(N) 



and 



(N) 



C/fc(cos((?!)^ - 0^)) 



for < A; < — 1 and < r < — 1 . The matrix m['V is the coefficient matrix 



of (3.5 1 when 7](t) = 1. We note that these are symmetric matrices. 
Theorem 3.3. For Q < k,r < N - \, 

(a) the matrix M^^-* is nonnegative definite with all eigenvalues in [0,1]; 

(b) the matrix Mjf^^ is positive definite if and only if k + r < TV; 

(c) If k + r > N, then zero is an eigenvalue of m\!^J which has multiplicity equal 
tok + r + l- N. 
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Proof. We start with an observation. Let fc = N — I ~ 2. Since (p^ = 2tti'/N, it 
follows readily that sin(fc + 1)(0^ — (p^) = — sin(Z + 1)(0^ — (p^). Hence, ii fi v 
then = —b''^^^, whereas b'^^, = k + l = N— {l + l) = N — bl,^. Consequently, 
we see that 

N 



(3.6) <i-2.,. = Ir 



NIr - <^ 



1 r(^) 



forO<iV-;-2<iV-lorO<;<Af-2. Thus, we only need to consider B^^J . 
Let us define column vectors coSj and siuj by 

cosj = (cos j(/>^)po and sin^- = (sinj>^)^^J„ j > 1, 

and let 1 = (1, . . . , 1) also as a column vector. It is well known that Un{t) can be 
expressed as 

C/2to(cos e) =2 cos 2me + 2 cos(2m - 2)6' + . . . + 2 cos 261 + 1 
U2m+i (cos 9) =2 cos(2m +1)^ + 2 cos(2m - 1)^ + . . . + 2 cos 6*. 

Using the fact that cosj(0;i — ^v) = cos jcp^ cos j (pi, + sinj^^ sinj(^,^, we can then 
write the matrix as 

^2^,r = 1 • 1^ + C0S2 • cos^ + sin2 • siug + . . . + cos2m • cos^„ + sin2m ' siu^^ 

= ^2171X2^, 

where ■= (1, COS2, sin2, . . . , cos2to) sin2m) denotes the matrix that has l,cos2, 
sin2, . . . , cos2„i, sin2m as its column vectors. In the case of A; = 2m + 1, we have 

-^2m+l,r ~ -^^m+lX-^m+li ^2rn+l '■= (cOSi , siui , COS3 , . . . , COS2„i+i , Sin2m+l) • 

Considering the quadratic form b'^^Jc, if necessary, this shows that the matrix 
Bff^J, hence N~^Blf^J = Ir — ^'kr\ is nonnegative definite. Consequently, we see 

that the eigenvalues of M^j^J are all bounded by 1. Furthermore, the matrix is 

of the size r x (fc + 1) so that its rank is at most min{A: + 1, r}. Consequently, if 
^2mC = for a vector c e then the trigonometric function 

T2m{i) '■= ci + C2 COS 2t + C2 sin 2t + . . . + C2m COS 2mt + C2m+i sin 2mt 

vanishes on the points t = (p,^ for < v < r — 1. If r > 2m + 1 = A: + 1, then the 
trigonometric polynomial T2m of degree k vanishes on at least 2m + 1 points, which 
implies that T2m{t) = 0, so that c = 0. It is easy to see that the same also holds for 
k = 2m + 1. Consequently, the columns of are linearly independent if r > k+l. 
If r < k+l, then we consider the rxr matrix, 1^, formed by the first r-th columns 
of Xk- Considering Y^c = as above, we see that Yfe has full rank. Consequently, 
rank(Xfe) > rank(yfe) > r. Thus, we have proved that rank(Xfe) = min{/c + l,r}. 

If A; + 1 > r then, for c e W, c^B^^Jc = {c^Xkf = so that c = as 
rank(Xfc) = r. This shows that b\^^ is positive definite, hence invertible. Whereas 
if fc + 1 < r, then the rank of b'j^^ satisfies 

rank(S^^^) > rank(Xfe) + rank(Xfe) - (fc + 1) = fc + 1, 

which shows that rank(i?|,^^) = fc + 1. Hence, B^^ is singular in this case. Con- 
sequently we have proved that b\^^ is positive definite if and only if fc + 1 > r. 
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Hence, by |3.6[ the matrix M^j^J is invertible if and only if N - 
is equivalent tofc + r + l<A^. 



fc — 2 + 1 > r, which 



Furthermore, if k 



1 < r, then the kernel of the matrix S^^'' has dimension 
r — (fc + 1) . It follows that zero is an r — (fc + 1) fold eigenvalue of the matrix. Again 



by (3.6), this is equivalent to that is a fc 



1 — fold eigenvalue of 



(N) 
k.r 



□ 



Since we need to solve (3.5 1 for all k = 0,1, . . . , N — 1, the above result shows 
that the method will not work with ri(t) = 1 for any r > 1. The role that rj plays 
then becomes essential. 



Let us define by A^'^J the coefficient matrix of the system ( |3.5[ ), 



Ir 



sin(fc + - (f)^) 



where Ir is the identity matrix of r x r. This is also a symmetric matrix. 
Theorem 3.4. For < k,r < N - I, 

(a) if k + r < N, then the matrix A^'f^J is positive definite with all eigenvalues in 
(0,1]; 

(b) ifk + r > N, then the matrix A^j^J is positive definite if and only if t < 1 —r/N. 
Proof. Let us denote the eigenvalues of a matrix A by /ij [A] . By the definition, it 



is easy to see that ^j{Ir — A^^^) — ri{jj)fj,j{Ir — Af^"''), which imphes that 



(3.7) 



1 



^(#)+^(#v.(*C) 



If fc + r < A^, then fij{MjfJ) > for fc + r < A^ by the theorem, and (|3J| implies 
that 

M.(47) > ^(^)/^.(^C^) > ^(1 - §)t^M7) > 

since 77 is non-increasing. Thus, for k + r < N , the matrix A^f^J is positive definite. 



\ — N. By (3.7), A^^^ has 1 — ^ '"^^ eigenvalue of 



If fc + r > A, then Af^"^ is nonnegative definite and has zero as an eigenvalue 
of multiplicity k 4 

multiplicity fc + r + 1 — A^, and a']^^ is positive definite if and only if 1 — 7](-^) > 0. 
Since fc + r > A^, we have fc = A^ — r. A — r + l,...,A— 1. The assumption 
T < 1 — r/N implies then that > t for fc + r > Aand, consequently, 1 — ri{j^) > 
as 77 is strictly decreasing on [r, 1]. On the other hand, if r = 1 — r/N, then 



V{^) - Vir) = 1, so that Ai!^J 
and, hence, is singular. 



has at least one zero eigenvalue when k > N — 



□ 



As a consequence of this theorem, the matrices A^jf^J are all positive definite, 
hence invertible, if r < (1 — t)N, where t is the cut-off point in rj. Thus, the 
condition r < (1 — t)N becomes a necessary condition for the algorithm to work. 
The reason that it is not sufficient lies in the numerical analysis. Theoretically, this 
condition is sufficient for A^f!^J to be invertible, but these matrices can be severely 
ill-conditioned which render the algorithm useless. For a positive definite matrix, 
the conditional number can be defined as the ratio of its largest eigenvalue over its 
smallest eigenvalue; that is, if A is a r x r positive definite matrix with eigenvalues 
^io, . . .,Hr-i, then 

cond(A) := max u,/ min u,-. 

0<j<r-l 0<j"<r-l 
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By (3.7) and the proof of the last theorem, if k+r > N, then the smallest eigenvalue 



of is 1 - 77(1 -r/N), which can be very small when r is close to 1 — r/N, as rj 
is strictly decreasing on [r, 1]. Thus, it is necessary to take r away from 1 — t /N , 
or, in other words, choose r < 1 — r/N + e for some e > 0, to prevent the matrices 
A^j^^ become too ill-conditioned. On the other hand, when k < tN, we have 

^fcr^ ~ ^'^kV ™^ ^^'^ matrices Mjf^J can be severely ill-conditioned. Thus, we 
often have to choose r fairly small. 

The eigenvalues of a related matrix, C$, were studied by Slepian in |16j . where 



sin2(^ — v)^ 



When < <i> < tt/2, the eigenvalues of C$ are all between (0, 1) and the asymp- 
totic of the largest eigenvalue fiQ is given in |16j , which shows that 1 — /ig can be 
exponentially decay as r ^ 00 (for precise statement, see [T6j p. 1387] with the 
notation Afe(r, <!>)). If $ = (fc + l)7r/N, then we see that 

_ fc+lsin(A:-|-l)((/)^-0^) 



N 



which is similar to our 5^ j,. For fixed r, k and N sufficiently large, the matrix C$ 
with <f> = (fc + 1)/A^ can be regarded as a close approximation to B^J , so that 
the eigenvalues of C$ gives some indication to the eigenvalues of and hence, 

those of . However, a small perturbation in the entries of the matrix may 

lead to a large change in the eigenvalues; thus, it is of interesting to understand the 
eigenvalues of M^^^ itself. 

It should be mentioned that the matrix C$ and its eigenvalues are instrumental in 
deriving the singular values of the Radon transform ([8,i9j) as well as in completing 
data using singular value decomposition for the limited angle problem. 



3.3. Algorithms for limited angle problem. We now consider the limited angle 
problem for which the given data set is (3.1 1 and we assume that N is even. With 



simple modification, the method will work with odd N as well. 



Recall that for N being even, we use (2.10) instead of (2.5), so that we replace 



N in the systems of linear equations in (3.5) by N/2 and the coefficient matrices 
of these systems are non-singluar, according to Theorem 3.4 if r < 1 — 2r/N 
provided Nd = N/2. Below we sum up the algorithm for limited angle problem and 
we assume Nd = N/2. 

Algorithm 3.5. (Algorithm for limited angle problem) Given Radon data 
{91^, k '■ r < II < N/2 — 1, < fc < N/2 — 1}, where N is an even integer. 

Setp 1. For li^r,..., N/2 - \, compute for k = 0,1, ... , N/2 by FFT 



N/2-1 



j=o 



(2j + 1)^ 

N 
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Step 2. For a given r choose r so that r < 1 — 2r/N and choose an rj. For 
fc = 0, 1, . . . , N/2 — 1 solve linear system of equations 

r-l N/2-1 

(3.8) Afe,^-^a^'lAfc,, = ^I'-^^k,., 0</i<r-l, 

for X^^kj < /X < r — 1, where 

S'tep 5. Augmenting Xk,u computed in Step 1 and Step 2 to obtain a full set 

Aw := {\k,^ : < J/ < iV/2 - 1, < fc < N/2 - 1} 
and applying OPED Algorithm\2.4\ on A^r to reconstruct the image. 



The output of the second step of the algorithm gives approximation for the 
missing data Ao,fc, . . . , Ar-i.fe for fc = 0, 1, . . . , N/2 — 1. Notice that the algorithm 
does not complete the data set itself, what it completes is the set of sine transform 
sXk,fi of the data. 

We now turn to the problem of how to choose rj. Let hk{t) be a polynomial of 
degree 2k + 1 such that hk{0) = 1, h\;!\o) = for 1 < j < fc, and h['\l) = for 
< j < k. Such a polynomial is given explicitly by 

hk{t)^{i-tr^Y.(''~.'V- 



J 



For a fixed k we then define rj^t) by 

fl, 0<t<T, 

(3.9) rj{t) -.^ [{E^) , T<t<l 

[o, t>l. 

Then rj E C'^(M) and it satisfies the desired property. The function rj curtails the 



values of high degree proj;. / in the expansion (2.4 1. Note that r] does not have to 
be zero at i = 1. In fact, we can choose r/ so that it is smooth on [0, 1], 77(1) = 1 
for < t < T and ri{t) decreasing to 77(1) — P > on [t, 1]. For example, here is 
such a function in C^, 

hkAt) (/3-l)(3i2-2i3) + l, 
which when used in (3.9 1 gives a function in so that 7^(1) — (3. 



Naturally then we face the problem of how to choose r and (3. As the discussion 
at the end of the previous subsection shows, we should choose r reasonably small 
to avoid the ill-conditioning of the matrices. The condition r < 1 — however, 
is only a necessary condition; we need, in practice, r substantially smaller. There 
is, however, a balance, as the algorithm preserves polynomials up to degree riV^. 
Small r means lower degree of polynomial preservation and less accuracy in recon- 
struction. This is where /3 comes into the picture. If j3 is large, say j3 = 0.95, then 
rj will decreasing slowly down from 1 to 0.95, and we will have almost polynomial 
preserving property. The experiments have shown that larger (3 may lead to worse 
condition numbers of the matrices, but the increasing is not drastic. On the other 
hand, the condition numbers increases drastically as r increases. 
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For a fixed N we can compute the condition numbers of Aj^ J numerically. We 
give an example. Notice that when r is fixed, the available data {gk.u '■ r < v < 
N/2 — 1, < fc < N/2 — 1} is over an arc of tt — 2Tir/N radiant or the missing data 
is over 

a := 2'Kr/N = (360r/iV)°. 

In other words, the given data is limited with angles over an arc of 180 — a degree 
and the missing data is over a degree. 

Let us take for example N = 502, which means the full data consists of 251 views 
of equally spaced angles over [0, tt] and 251 rays per view. For the incomplete data, 
if r = 21, then the available data is limited to an arc of 165°, a 15° difference from 
the full data. If r — 42, then the data is limited to an arc of 150°, a 30° difference 
from the full data. In Table 1, the the maximum of the condition numbers for our 
matrices, rounded to nearest integers, are given for different values of r and (3 in 
the cases of r = 21 and ?■ = 42. 



Table 1 . Maximum of condition numbers 



r = 21 









r= 42 








T 




max 




T 




max 


0.0 


0.5 


44 




0.0 


0.5 


135 


0.0 


0.9 


160 




0.0 


0.9 


503 


0.1 


0.5 


293 




0.1 


0.5 


60295 


0.1 


0.9 


716 




0.1 


0.9 


68296 


0.2 


0.5 


48900 




0.2 


0.5 


3.66715 X 10^° 


0.2 


0.9 


48928 




0.2 


0.9 


3.66715 X 10^° 



For example, in the case of r = 21, r = 0.0 and (3 = 0.9, the maximum of the 
condition number is merely 160. The maximum is very large in the case of r = 42 
and r = 0.2, showing that the matrix A^^^ is severely ill-conditioned for some k in 
this case. Furthermore, the maximum of the condition numbers appears to increase 
drastically as r increases as well as r increases. Another interesting fact is that the 
dependence on (3 appears to be insignificant for larger r and larger r. In the Figure 
[l] the distribution of the condition numbers in the case of r = 42, t = and r = 0.2 
is plotted, which shows that not all matrices among a'^^^ become ill-conditioned. 



v. 



Figure 1. Condition numbers for r = 



42. Left: r = 0. Right: r 0.2. 
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An interesting fact is that the conditional numbers in the case of r = remain 
reasonably in check even when r is large, as seen in the following table, where we 
choose P = 0.9 to compensate t = 0. 

Table 2. Maximum of condition numbers for t = and (3 = 0.9 



r 


21 


42 


63 


83 


126 


max 


160 


503 


1037 


1757 


4084 


limited angle 


165° 


150° 


135° 


120° 


90° 



In the case of r = 126, the given data is distributed over an arc of 90°, which 
means that half of the full data. In this case, the maximum of the condition 
number is 4084 for (3 = 0.9, which is still not too large. However, r = means 
that the algorithm no longer preserves polynomials and this is the case that should 
be avoided. Still, by choosing (3 large so that the result of the sampling on the 
coefficients is not too far away from polynomial preservation, the case t — can 
be used to reconstruct of images as our numerical tests have shwon. In general, 
however, we should work with positive r whenever we can. This is supported by 
the numerical experiments discussed in the next section. 

4. Numerical Experiments and Discussions 

We have applied the algorithm in the previous section on several examples, which 
are presented and discussed below. Recall that our data of limited angle consists 
of {gv,k : r < V < N/2 — 1, < fc < iV/2 — 1}, where N is an even integer, and the 
angle within which the data is distributed is, for a given r, 180° — (360r/A^)°. 

4.1. Shepp-Logan phantom. For our first numerical example, we use the classi- 
cal head phantom of Shepp-Logan [T7] . This phantom is shown in Figure [2] 




Figure 2. Reconstruction based on full data 



The left figure is the original phantom. The right figure is the reconstruction 
by OPED based on the full data with N — 502, which means 251 views with 
angles equally distributed over [0,7r] and 251 rays per view, and the size of the 
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reconstruction is 256 x 256 pixels. Reconstruction based on the full data has been 
discussed in [H [Ml HI] , we will not give further details here as our purpose is to 
demonstrate the feasibility of our method on the limited angle problem. 

For the reconstruction on the limited angle data, we choose the same set-up, 
with 201 angles over [0,7r] and 201 equally spaced parallel rays in each view. 

In our first example, r = 21, which amounts to data limited in an angle of about 
165°; in other words, views from about 15° angle are missing. The reconstruction 
by our algorithm is given in Figure [3] in which (3 — 0.9 and r = for the left figure 
and 0.2 for the right figure. 




Figure 3. Reconstruction with r = 21. Left: r = Right: t = 0.2 



The left image is reconstructed with r = and /? = 0.9; it is a fairly accurate 
reconstruction, although there are noticeable artifacts in the direction of missing 
views and a bit distortion around two spots on the edges. The right image is 
reconstructed with r = 0.2 and f3 — 0.9; it shows clearly artifacts of ripples, but 
the image appears to be sharper and has less distortion than the one in the left 
otherwise. In the case of t = 0, the maximum of the condition numbers of the 
matrices A^f^J is 160, so that the matrices are rather well conditioned. In the case 
of T = 0.2, the maximum of the conditions numbers is 48928, which may have 
contributed to the ripples in the image. 

The condition that guarantees the non-singularity of the matrices in this case is 
T < 1 — 42/502 « 0.916335, whereas our computation of eigenvalues shows that r 
has to be much smaller in order that the matrices are well conditioned. For our 
other examples, we will mostly take r = 0. The choice of /? = 0.9 means that our 
sampling of coefficients follows a curve that decreases from 1 to 0.9, a decline that 
is rather mild, which leads to reasonable reconstruction image. 

In our next example, we consider the case r = 42, which means that the data 
is limited to views with angles distributed over an arc of 150°. The reconstruction 
with r = and /3 = 0.9 is given in Figure |4] 

In this image, artifacts and distortion are clearly visible and most prominent at two 
points on the edges of the images. The maximum of the conditional numbers in 
this case is merely 503, so that the matrices are in fact fairly well conditioned. This 
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Figure 4. Reconstruction when r = 42. 



suggests that the distortion is hkely caused by the choice of t = 0, which means 
that no polynomial preservation is kept. 

4.2. Data with noise. The limited angle problem is well known to be ill-posed. 
Below we present our reconstruction with noise data. We use again the Shepp-Logan 
head phantom but add noise in the data, which is Gaussian normally distributed 
with zero mean and a standard deviation 0.03. The noise is about 2% in the data. 
For limited angle, we choose r — 21 and 42, respectively, which correspond to data 
limited over an arc of 165° and 150°, respectively. The reconstructed images by 
our algorithm with r = and (3 — 0.9 are given in Figure [5] 



Figure 5. Noise data. Left: r = 21. Right: r = 42 

These reconstruction should be compared with the left image in Figure 3 and the 
image in Figure 4, respectively, which are the reconstructed images based on the 
same limited angle data but without noise. These images indicate that our method 
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is relatively stable, in the sense that the reconstructed images are not distorted 
much by the noise. 

4.3. Discussion. The theoretic study and the numerical experiments point out 
that the proposed algorithm depends critically on the choice of r. The matrices 
remain relatively well conditioned for t — even when r is large, but the case t = 
introduces distortion in the images, in addition to the artifacts. The reconstruction 
with T > appears to lead to less distortion in the images. However, the maximum 
of the condition numbers appears to grow exponentially with r for r > and 
it increases drastically still for larger r. The ill-postcdness of the matrices likely 
reflects the ill-posed nature of the limited angle problem. It is likely that solving 
the linear systems with pre-conditioning algorithms may improve the reconstructed 
images. This is, however, beyond the scope of the present paper. 

5. Conclusion 

A method for reconstruction images in the limited angle problem is presented 

and a theoretic study is carried out. The ill-posed nature of the problem shows up, 
when r is not zero, in the ill-condition of the linear systems of equations that we 
need to solve. Numerical tests have demonstrated the feasibility of the method. 

In order to fully understand the proposed method, further numerical study needs 
to be carried out. One interesting question is how much of the artifacts and the 
distortions are due to the ill-conditioning of the matrices when r is not too small. 
The theoretic study indicates that the algorithm should be applied with t relatively 
large if the severely ill-conditioned systems can be solved. On the other hand, as 
the limited angle problem is intrinsically ill-posed, there will have to be distortion 
of images when the angle is small. 
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